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Abstract 

Recent work by McClarren & Hauck [29] suggests that the filtered spherical harmonics method represents 
an efficient, robust, and accurate method for radiation transport, at least in the two-dimensional (2D) case. 
We extend their work to the three-dimensional (3D) case and find that all of the advantages of the filtering 
approach identified in 2D are present also in the 3D case. We reformulate the filter operation in a way 
that is independent of the timestep and of the spatial discretization. We also explore different second- and 
fourth-order filters and find that the second-order ones yield significantly better results. Overall, our findings 
suggest that the filtered spherical harmonics approach represents a very promising method for 3D radiation 
transport calculations. 

Keywords: Radiation transport, Pat— method, spherical harmonics, asymptotic diffusion limit, 
discontinuous Galerkin 



1. Introduction 

Many phenomena in the universe involve the transport of radiation and need to be modeled with 
radiation-transport techniques that are as accurate as possible to maximize the match with observations. 
Examples are nova and supernova explosions, gamma-ray bursts, star or planet formation, luminous blue 
variable outbursts, stellar winds, etc. In these examples, radiation plays a major role in exchanging energy 
and/or momentum between different parts of the system. In most of these cases, the radiation is composed 
of photons, but the radiation can also be composed of neutrinos in studies of core-collapse supernova explo- 
sion mechanisms (see, e.g.^ [3 US]) or in modeling the torus orbiting the black hole produced in the merger 
of neutron stars (see, e.^., [42]). 

One of the main difficulties in solving the radiation-transport equation arises from the multidimensional- 
ity of the problem. Radiation is described not only by the position of the radiation carriers, but also by their 
direction of propagation and energy (or, equivalently, by their momentum), making the problem (6 + 1)- 
dimensional in the most general case (3 dimensions for the spatial coordinates, 2 for the angular direction, 1 
for the energy, and 1 for the time). Another source of complication stems from the fact that many problems 
consist of regions of strongly-varying optical depth. For example, many astrophysical systems contain a 
central source of radiation, where the optical depth is high, surrounded by more transparent outer regions. 
Due to high opacity near the center, radiation migrates out of that region mostly via diffusion, while in the 
outer parts, it streams to infinity more freely without much interaction with matter. Correspondingly, in 
the limit of infinite optical depth, the transport equation acquires a parabolic character, while it maintains 
a hyperbolic one in all other regimes (see, e.^., [33 ). Radiation-transport approaches must therefore handle 
accurately not only these two distinct regimes, but also the transition between the two, which is generally 
the most difficult to treat. These features make the solution of the transport equation both complex and 
computationally expensive. 

There are two commonly used ways of simplifying the solution of the transport equation. One approach 
consists in reducing the number of degrees of freedom by assuming spherical or axial symmetry. While this is 
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a good simplification for some problems, there are many other situations in which the system does not posses 
any spatial symmetries and hence the transport equations need to be solved in three spatial dimensions. 
Another way of simplifying the problem is via approximating the form of the transport equation (this 
is equivalent to reducing dimensionality in the momentum space). One of the simplest examples is the 
diffusion approximation, where one approximates the transport equation with the diffusion equation {e.g., 
[39 ). This makes the equation far simpler and computationally less expensive to solve {e.g., ES])- However, 
there is a price to pay for this: Although the diffusion approximation is accurate at high optical depth, 
it leads to incorrect results in low optical depth, i.e. "transparent", regions. This can be improved by 
using flux limiters {e.g., [14 ), but the diffusion approximation still cannot correctly capture the anisotropy 
of radiation commonly encountered in those regions (see, e.g., [38]). Moreover, the treatment of semi- 
transparent regimes is somewhat artificial because the radiation fluxes in those regions are usually calculated 
via interpolation between the two fluxes calculated assuming free- streaming and diffusive transport. There 
are several more accurate ways of approximating the transport {e.g., two- moment schemes with analytic 
closures, etc.; e.g., [8 ), but the solution of the full (6+l)-dimensional transport equation is often not only 
desirable, but actually necessary for an accurate modeling. 

One of the most commonly used methods for solving the transport equation in the multidimensional 
(both in space and momentum) case is the discrete-ordinate {Sn) method, which solves the transport 
equation along several directions in each spatial zone [Ml [38l HH [20] . Unfortunately, this method has several 
drawbacks. Most importantly, it suffers from "ray effects" in multidimensional calculations [34]. Due to the 
discrete nature of the angular representation, in fact, radiation cannot reach regions between these discrete 
directions, leading to large spatial oscillations in the transport variables. In addition, time-implicit Sn 
methods require very complex solutions and parallelization procedures [4]. 

Monte-Carlo methods {e.g., [T8l|19l[2J) are often regarded as the most accurate method for radiation 
transport, but they are also not without drawbacks: Monte-Carlo solutions exhibit statistical noise due to 
the finite sampling of the phase space. Since this noise decreases only as A^~^/^, where N is the number of 
Monte Carlo particles, it can take many particles to produce a sufficiently smooth solution, making large 
simulations computationally very expensive. 

Another approach to radiation transport is the spherical harmonics, or P/v, method. This method is based 
on an expansion of the radiation intensity (or distribution function) in angles using spherical harmonics basis 
functions. This results in a hyperbolic system of partial differential equations for the expansion coefficients, 
which represent angular moments in this basis. The spherical harmonics method has several interesting 
characteristics. Due to hyperbolicity, the P/v equations approximate radiation as a series of waves with 
velocity bounded by the speed of light [30] . This restriction is consistent with the transport equation, 
in contrast to diffusion methods where radiation propagates at infinite velocity. Moreover, the spherical 
harmonics expansion exhibits formal spectral convergence to the true solution. Such an expansion also 
preserves the rotational invariance of the solution, unlike Sn methods, where the absence of such invariance 
leads to the ray effects mentioned above Another advantage of the P/v method is that it generally uses less 
memory to model a given angular distribution at a given accuracy compared to, e.g., the Sn method. A Pn 
truncation is roughly equivalent to a Sn-\-i solution, but, while the former has (A^ + 1)^ degrees of freedom, 
the latter has 2{N + 1)^, thus it requires roughly twice as much memory. Given that memory requirements 
represent one of the main difficulties in 3D radiation-transport calculations, a factor two in memory saving 
represents a significant advantage. 

However, spherical harmonics methods also have some negative aspects. Most importantly, in transparent 
regions, the solution to the P/v hyperbolic system exhibits non-physical oscillations. These oscillations 
are related to the so-called Gibbs phenomenon that occurs when non-smooth functions are approximated 
with smooth basis functions [7]. The worst consequence of such oscillations is that they can cause the 
radiation intensity to become negative, which may lead to negative matter temperatures when the radiation 



■"^Here, rotational invariance means that the operators of angular discretization and of rotation in space are commutative. 
In other words, the result of any rotation and then of a spherical harmonics angular discretization is the same as that of an 
angular discretization and then of a rotation 7 . In the method, this is true only for those rotations which map the angular 
grid onto itself. 
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transport is coupled to the matter energy equation [30l |35] . There have been several attempts to address 
this problem ^3 ^3 |9l [29l |36]. One of the most efficient, robust, and accurate approaches is the one by 
McClarren & Hauck [29], in which filters are proposed to remove oscillations of the radiation intensity. By 
filtering out the oscillations, McClarren & Hauck [29 were able to avoid negative solutions while maintaining 
high angular accuracy. Their filtering approach also has the advantage of being easily extendable to high- 
order P/v expansions, producing formal convergence to the transport solution and preserving the equilibrium 
diffusion limit. 

Although the results of [29 strongly suggest the idea of applying a filter to spherical harmonics expansions 
is an efficient, robust, and accurate way of doing P/v radiation transport, there remain several open questions. 
For example, [29 considered only one type of filter, the so-called spherical-spline expansion [7]. There are 
other types of filters that have some interesting theoretical and numerical properties [7 , which might be a 
more optimal choice for application to Pn transport. Moreover, McClarren & Hauck [29 presented results 
only for the 2D case, leaving open the question of how well their filtering approach performs in 3D. Also, 
this filtering scheme, as realized by [29 , and as we will discuss in detail later in this paper, does not have 
a clear continuum limit as the spatial resolution and/or timestep approaches zero, which implies that the 
solution cannot be studied for spatial convergence. 

In this paper, we reconsider the P/v scheme and extend the work by McClarren & Hauck [29 in a number 
of ways. Firstly, we reformulate the filtering procedure in such a way that it acquires a clear continuum 
limit. Secondly, we investigate a wider range of filter types and strength parameters. Finally, we perform 
calculations both in 2D and 3D. For this we have we have developed the new radiation transport code 
Charon, whose ultimate goal is that of performing 3D, general relativistic multi-energy, multi- angle and 
velocity-dependent radiation-transport calculations. In this paper, we present the ffist step towards this 
goal. 

Charon uses the semi-implicit scheme of McClarren et al. [28 for time integration. In this scheme, the 
streaming parts are treated explicitly with a second-order Runge-Kutta method, while the matter-coupling 
terms are treated implicitly because of the stiffness introduced by the coupling. Also, the implicit system 
of [28 is local to each spatial element, and for linearized matter-coupling terms this implicit integration 
scheme becomes trivial [28 . 

Although the timestep in the semi-implicit approach is limited by the Courant-Friedrichs-Lewy (CFL) 
condition based on the speed of light, this is not a serious drawback for the kind of applications that we have 
in mind. In radiation-hydrodynamics simulations, the timestep size would still be limited by the dependence 
of the matter properties (e.^., opacities and emissivities) on the matter temperature (and electron fraction, 
if we are dealing with neutrinos with conserved lepton number). Moreover, the semi-implicit approach is 
relatively easy to parallelize {e.g., via domain decomposition) and has lower communication requirements 
compared to fully implicit schemes. This should translate to significant advantages in massively-parallel 
large-scale radiation-hydrodynamics simulations. Moreover, in many relativistic systems, the sound speed 
of the fluid is comparable to the speed of light, and thus the radiation and hydrodynamics timescales 
are comparable, reducing the extra cost in an explicit treatment of the radiation streaming. The spatial 
discretization in Charon is based on the asymptotic-preserving linear discontinuous Galerkin (DC) scheme 

This is the ffist in a series of three publications, with the present one dealing with static matter con- 
figurations and focusing exclusively on the transport methodology. For this reason, we do not consider the 
coupling with hydrodynamical equations. As a result, when an absorption occurs, the particles are sim- 
ply removed from the system and do not increase the matter temperature. The second paper will include 
velocity-dependence and coupling to hydrodynamics, whereas in the subsequent publication we will present 
a fully general-relativistic algorithm. 

The paper is organized as follows. In Section [2j we introduce the concept of radiation distribution 
function and the relativistic Boltzmann equation. In Section [3j we describe the numerical schemes used in 



our code for frequency (Section 3.1), angular (Section |3.2[ ), spatial (Section |3.4[ ), and time discretization 
(Section |3.5| ). In Section [4j we present numerical tests of these schemes. Finally, we summarize our results 
and provide our conclusions in Section [sj We use a spacetime signature (—,+,+,+), with Greek indices 
running from to 3 and the Latin indices from 1 to 3. We also employ the standard convention for the 
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summation over repeated indices. 



2. The relativist ic Boltzmann equation 

2.1. The distribution function for radiation 

Radiation is usually described in terms of the specific radiation intensity, /, defined such that 

dS' = I cosOdAdudQdt, (1) 

represents the energy of radiation in frequency range du centered around traveling in direction Q confined 
to a solid angle element dQ, which crosses, within time interval dt, an area dA oriented such that is the 
angle between the normal to the surface dA and direction Q {e.g., [39 ). In the case of neutrino transport 
and/or in the relativistic case, it is more convenient to work directly with the distribution function, F, which 
gives the density of radiation carriers on a given point in phase space. The reason for this is that (1) the 
distribution function is a Lorentz-invariant quantity [33 , and (2), as we discuss in more detail later, the 
distribution function allows us to compute the number density and the energy density of the radiation in a 
more natural way. 

In order to introduce the distribution function, we first discuss the notion of single-particle phase space 
in special relativity. In this picture, the particles are described in terms of their positions in spacetime, 
x^, and their momentum four- vector, p^, as measured in a fiducial inertial frame. Using the normalization 
condition for timelike vectors, the four- momentum has only three independent components, which can be 
expressed in terms of three spatial components, p\ or in terms of radiation frequency, z^, and two angles 
(6>, (/)) that describe the propagation direction}^ 

= ^ (1, cos^sin6>, sin sin ^, cosO) . (2) 
c 

Since we wish to define the distribution function in terms of p% or, equivalently, in terms of u^cj) and 0, 
we need to construct a Lorentz-invariant volume element, dll, over the manifold of the allowed momentum 
four- vectors. This is accomplished with the choice [15] 

dn = ^i^^^^ = ^d.di^. (3) 

-Po 

The distribution function is then defined in such a way that the quantity 

dN = F pH^ d^x dU = F d^x dv dVt , (4) 

is the total number density of radiation carriers in a spatial volume element d^x and phase-space volume 
element dll with trajectories traversing a t = const hypersurface with normal t = Q^j^ 

Since d^x = dA cos 6 dt and the energy per particle is given by hv, we have d(o = hu dN. Using this and 
comparing equation ([T]) with Q, we obtain: 

I=^F. (5) 
Note that since dA^ is a scalar, F is also a scalar quantity]^ 

^Assuming the radiation carriers to be massless, their energy and frequency are simply related as £ = hu. 
^Note that d'^x is not a Lorentz-invariant quantity, while p^t/j.d^x is one. 

^Note that this relation slightly differs from the one frequently encountered in the neutrino-transport literature {e.g., 11): 

/ p p 

where g is the statistical weight of the particles (g = 1 for massless neutrinos, g = 2 for photons) and e is the particle energy. 
This is due to three reasons: First, our specific intensity given by Eq. is defined in terms of energy per frequency interval, in 
contrast to energy per energy interval, as used in the neutrino-transport literature. Moreover, our distribution function already 
contains the factor g, as can be seen from its relation to the total number density of radiation carriers given by Eq. (Q. Finally 
we use the Lorentz-invariant volume element given by Eq. jsj instead of d^p. 
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2.1.1. The relativistic Boltzmann equation 

The special-relativistic Boltzmann equation can be written as [33] 



where C is the cohisional term describing the interaction of radiation with matter, while the left-hand-side 
of the equation describes the propagation of radiation. In order to compute C, we express it in terms of the 
absorption, emission and scattering coefficients. To do that, we start by considering the evolution equation 
for the intensity of radiation [39], 



c dt ^ dx^ ^ 



where 77 represents the radiative emissivity of the matter, n is the total extinction coefficient and combines the 
absorption and scattering coefficient^ and and K is the scattering kernel, expressing 

the probability of scattering from a given angle and frequency over to another angle and frequency [39 . 

Using equations ([2|, ([5|, and ([7|), it is easy to obtain an equation for F in terms of the above extinction 
coefficients: 

^'£=5^-'^^'^^+^ / (7)W'^Piw/^d.'di^'. (8) 

Notice that, since F is a scalar, so is C[F], thus we find the classical result that rf/v'^^ and vn are invariant 
quantities [33 . 

Our scheme is in principle able to handle any type of scattering kernel, but for simplicity, here we will 
only consider the case of elastic scattering, e.g.^ scattering in which the radiation energy does not change. 
In this case, the scattering kernel can be expresses as pT] . 

K{v',n' ^v,n) = ^[l ^ (Jafi ■ n']8{v - v') , (9) 
47r 

where the scattering anisotropy is modeled using only one coefficient Ga- 



3. Description of the scheme 

In general, the distribution function, F, is a function of 7 variables: the time and spatial coordinates, 
x^, the frequency v and the angles of propagation ip and 0. These variables are usually defined either in 
the Eulerian (inertial) frame or in the comoving frame (ie., a set of frames, each of which has a velocity 
that instantaneously equals that of the matter element, e.g.^ [33l|24]). In the case of static matter, as the 
one considered in this paper, these two frames are identical. In the scheme implemented by the Charon 
code, the distribution function is expanded in the spatial coordinates using the linear DG basis and in the 
angular variables using spherical harmonics. The frequency is treated using the multi-group approach. This 
yields a large system of ordinary differential equations that is then evolved in time using a semi-implicit 
time integrator. The details of the discretization are discussed in this Section. 

3.1. Frequency discretization 

We consider the case in which the distribution function has compact support in a frequency space 
given by the interval [0, i^max]- Although this is not strictly valid in the general case, radiation usually has 
negligible contribution above some cut-off frequency. Therefore, in many practical applications, one can 
choose i^max to be sufficiently large so that there is little radiation beyond this frequency. For simplicity 
of illustration, we introduce a uniform grid in this frequency space as Vn = nAu^ n = 0, 1, . . . , TVjy + 1, 



^Hereafter, the absorption or scattering extinction coefficients are defined as absorption or scattering opacitites or inverse 
mean-free paths. 



where = i^max/{J^u + 1) (an extension to a non-equidistant grid is conceptually trivial). The associated 
intervals [z/^^^n+i] are commonly called frequency or energy groups. Using these groups, we can construct 
an orthonormal basis {Xn}n=o 

/A fl/VK, if G [i^n,i^n+l] , r^^\s 2. / 3 

Xi(^) = S^ . . Vn= /I'^z/ dz/ -<). (10) 

10, otherwise, Jj^^ o 

We then expand a function / G i^^(0, z^max) on this basis as (for clarity we report the summation symbols 
in the expressions below) 

/iv.(^)=E/"^-(^)' f{u)h'u'du. (11) 

The truncated expansion, /at^, is then a first-order accurate (in L^— norm) approximation of /. 

S.2. Angular discretization 

As orthonormal basis on the unit 2-sphere, we consider the real spherical harmonics, (see 



Appendix A), whose orthonormality conditions are 



As a result, any function / G L'^{Si) can be expanded in spherical harmonics as 

N £ „ 

/^((^,^) = ^ ^ /-r,^((^,^), f{^,o)Y'^{^,o)dn, (12) 

where we have used the notation Y^^ to denote the complex conjugate of Y^m (which is equal to Y^m since 
we are working with real spherical harmonics). If / is a smooth function, will converge to / with spectral 
accuracy in the L^— norm [7]. 

3.3. The multi-group P/v scheme 

We consider the following ansatz for the expansion of the distribution function: 

F{x", „,^,e) = Y,T. E Xnii^) YeU^, e) , (13) 

n=0 ^=0 m=-£ 

where 

Fn^m(^")= / h^V^dv I d^lF{x'',V,^,e)Y,m{^,e)Xn{^)- 

Jo JSi 

To simplify the notation, we introduce the multi-index A = {n,^, m}, and the basis functions, 

^a(^, V^, 0) = Xn(^) Yim{^, 6) , 



SO that Eq. (13) becomes 

F(x", e, ^,e) = Y, F^{x^)^A{e, 9) = F^^ a • (14) 

A 

Note that the space spanned by our basis a} is a vector space. We adopt the usual convention of denoting 
vector components with an upper index and co-vector components with a lower index. Linear operators 



acting on vectors and co-vectors will have upper and lower indices associated with their decomposition on 
an appropriate tensor product combination of the canonical basis {^a} and its dual {^^}. 
Inserting Eq. (14) into Eq. (Ig]), we obtain: 



B +P 



C[F]. 



(15) 



Multiplying Eq. (15) by ^ and integrating with respect to dll, we then obtain 



dt ' ' ^ dx^ 

where we have again used Eq. ([2| and exploited the orthonormality of the basis. We have also defined 



(16) 



r B 



dn, 



and 



§^[F] = J C[F] ^^dn. 



(17) 



(18) 



The coefficients (17) can be computed exactly using a quadrature formula of high-enough order (see 



Appendix A). Since they are independent of position and time, they can be pre-computed and stored 
for later usage. 

The spectral decomposition of V^^b-, which determines the behavior of Eq. (16), is well known, see 
e.g. [10]. In particular, it has been shown that the eigenvalues are strictly bounded by the speed of light 
c. While this implies that there is no superluminal propagation of radiation, it leads to slower-than-light 
motion of radiation waves for finite N (the radiation velocity converges to the correct value with increasing 
N) [9]. This is particularly evident in low-order Pn free-streaming solutions. For instance, the maximum 
propagation speed for Pi is c/>/3. Filtering can also affect the propagation velocity of radiation [29 . In the 
multidimensional case, V^^b has also zero-speed modes that have to be treated carefully in Godunov-based 
schemes to avoid numerical instabilities [10 . 

In the general case, the source term (18) has to be computed at run-time, but for the particular case 
of a source term in the form ([o]), and assuming that the opacity coefficients are constant in each of the 
energy groups (as commonly done in multi-group schemes [39]), the source terms can be pre-computed up 
to constant factors. Under these assumptions, the source term becomes 



S^bF^., 



(19) 



wher43 
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l-'^B projector perpendicular to 



00, 



\ A — s:A 
^ B=^ B 



rA rnOO 
nOO^> L 



and A B is the anisotropy matrix: 



which can also be pre-computed. Note that we have denoted the opacity coefficients in the energy group 
[i^n, Vn+i] with the subscript n. 



7 





h-1/2 


k-\-3/2 










+ 




^ B * 


i B 





Xi-1/2 Xi Xi^i 

Ax 

r H 



Figure 1: The grid structure for the spatial hnear discontinuous Galerkin discretization. The white squares represent the cell 
centers, while the black dots show the cell interfaces. The dotted lines show the values of the Lagrangian basis li-i/2, h+2,/2- 
Finally the + and — show the interfaces where the inter-element fluxes, ^ see (|23|, are computed. 



3.4- Spatial discretization 

We discretize the system of Eqs. (16) in space using the asymptotic-preserving (AP) hnear DG scheme, 
see e.g. ^26|i28]. We recah that a scheme is AP if it reproduces a discretization of the diffusion hmit of the 
continuum transport equation in the hmit of smah mean-free-path. This is an important property since it 
guarantees that the diffusion of radiation has a correct rate even if the mean-free-path is smah compared to 
the spatial grid size. If a scheme is not asymptoticahy preserving, then the diffusion rate becomes unphysical 
when the mean-free-path is unresolved. 



For simplicity of notation, we consider a simplified ID version of Eq. (16): 



^ + ^ -^ = '^' 

where the multidimensional case will be discussed at the end of this Section. Furthermore, we employ a 
uniform numerical grid Xi = iSx^ while extension to a non- uniform grid is conceptually trivial. We then 
choose the points of the grid that are used to construct elements of width Ax = 2Sx, as shown in Fig. [l] 
This rather uncommon grid structure has been adopted to ease the integration of Charon with existing 
general-relativist ic hydrodynamic codes that use traditional finite- volume schemes. 

In the classical linear DG scheme, the degrees of freedom are usually identified with the one-sided limits 
of the solution at the points x^_i/2 and from the interior of the element, while in our case we evolve 

the cell-centered values at Xi and x^+i. The semi-discrete equations for the evolution of these cell-centered 
values can be easily obtained, since any function u in the finite-element space can be written as 

U{x) = Ui_i/2li-l/2{x) + ^i+3/2^i+3/2(^) , 

where 

h-i/2{x) = 1 , h+s/2{x) = , (20) 

— ^i-1/2 ^i+3/2 — ^i-1/2 



'Notice that the term in square brackets is the number of emitted particles. 
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so that 



and 



= ^^i-1/2 + ^^i+3/2 • 



(21) 



(22) 



The resulting numerical scheme is 



at 



where the flux terms in the element [xi_i/2, ^i+3/2] are given by 



3^+ 
2' 



(23) 



Here, 



is the average flux, while 



is the flux computed from the exact solution of the Riemann problem at Xi_i/2 with the left "L" and right 
"i?" states F"^L and F"^ji. The term is the flux across and is calculated analogously to . In 



the previous equations, we have decomposed V^^b as 



IC 



D 



C 



ID 



where IZ^^c are the matrices of the right eigenvectors, while h}^ d and C^^b are the eigenvalues and the 
left eigenvectors of V^'^b^ respectively. We have also introduced v > which is taken to be the first positive 
abscissa of the adopted Legendre quadrature and it is used to introduce extra numerical dissipation on the 
zero speed modes similarly to what is done in [10 . Notice again that the spectral decomposition of V^^Bi 
discussed in the previous Section, can also be pre-computed at the beginning of the calculations. 

Having described the numerical scheme for the ID problem, we can construct the multidimensional 



numerical scheme for Eq. (16) by repeating the same construction in every direction: 



dF^ 



dt 



Ax 



Ay 



i,j,k : 



(24) 



where the fluxes in the y and z direction, G and H, are computed analogously to the ones in the x direction. 

To avoid creation of false extrema in the numerical solution we use the slope limiting technique [17 . 
Among the different limiters that we have implemented are (1) the so-called "step limiter", which simply 
reduces the scheme to a first order discontinuous Galerkin scheme, (2) the "minmod" and (3) the asymptotic- 
preserving "minmod2" limiters [31 . The reason for using these particular limiters is that they have been 
well tested in the context of the transport equation (e.^., [31 ). Furthermore, the minmod2 limiter has been 
studied in detail in the context of linear DG methods, where it has been shown that it does not affect smooth 
solutions away from local extrema [16^ , thus yielding a scheme with a very small dissipation. 
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3.5. Time discretization 

For the time integration, we use the predictor-corrector method proposed by McClarren et al. [28]. In 
this approach, the streaming terms that model the transport of radiation are treated exphcitly, while the 
source terms responsible for interaction with matter are treated implicitly. The use of this particular time 
integrator is motivated by the fact that this yields a relatively inexpensive, stable and asymptotic-preserving 
scheme. As discussed in Section [l] the fact that the streaming terms are treated explicitly makes this 
scheme particularly easy to parallelize, while the associated CFL constraint is not particularly demanding 
for applications involving fluid moving at relativistic velocities and general-relativistic gravity. 



In order to simplify the notation, we rewrite Eq. (24) as 



"dT 



= e 



(25) 



where ^^[F] is a shorthand for the treatment of the spatial flux terms. For the time integration of Eq. J25b, 
we use the following two-step semi-implicit asymptotic-preserving scheme. Given the solution F^]^ at time 
/cAt, we first perform a predictor step 



F 



A 



fe+1/2 



At/2 



fe+1/2 



to obtain the solution at time {k + 1/2) At and then a corrector step: 



^ fe+1 - ^ k 

At 



= ^^[i^fe+1/2] + + S^sF^k+i 



to obtain the solution at time {k + l)At. At both stages, the absorption, emission and scattering terms 
are treated implicitly, while the streaming terms are treated explicitly. The explicit part of this scheme is 
second-order accurate in time, while its implicit part is first-order accurate [28] , 



3.6. Filtering 

Filtering is a common procedure to reduce the effects of the Gibbs phenomenon in spectral methods for 
numerical solution of partial differential equations [iTSl [40] . Filtered spherical harmonics expansions have 
been successfully used in, e.g., meteorology (see e.g. [7] and references therein) and the effects of filtering 
on the truncation error of a spectral expansion are now reasonably well understood [46]. However, the use 
of filters to mitigate (and, in most situations, remove) the occurrence of negative solutions in P/v schemes 
has only been proposed recently by McClarren and Hauck [29] . 

In their work, the authors propose to filter the spherical harmonic expansion of the solution after each 
timestep using a spherical-spline filter. Applying this suggestion to the spherical harmonic expansion of F 
we obtain (for clarity we report here the summation symbols) 



N 



e=o m=-l 



where 



cAt 1 



1 



(26) 



(27) 



Ax 7V2 (cr^x + 7V)2 ' 

and L is a characteristic length scale used to make a dimensionless, while at is chosen to be of the same 
order of magnitude as k,. 

The filtered spherical harmonics, FP/v, scheme has several interesting properties. Filtering has been 
found to be very effective and robust in removing numerical oscillations in P/v solutions, while preserving 
the rotational invariance of the scheme. Furthermore, for this particular choice of a, filtering turns off 
automatically in the limit N ^ 00, thus it does not spoil the convergence of the spherical harmonics 
expansion. 
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However, one important drawback, also remarked by McClarren and Hauck [31], is that the filtered Pn 
scheme does not have a clear continuum limit as Ax, At ^ 0. This is unfortunate because it implies that 
the filtered P/v scheme, FP/v, cannot be interpreted as a system of partial differential equations. This in 
turn means that the quality of the FP^ solution will depend on the spatial grid resolution in a way that is 
hard to predict. The ultimate and most important implication is that an FPn solution cannot be studied 
for spatial convergence. 

To solve this problem, we propose a modification/generalization of the FP^ scheme as follows. We 
introduce a strength parameter, 5 > 0, to be specified later, and construct the filtered expansion as 



N e r . ^ 



F'"'Yi^{ip,0), (28) 



where a{T]) is a fitter function of order p, that is, a function a gC^(M+;[0,1]) such thal^l 

a(0) = l, cr(^)(0)=0, for = l,2,...p- 1. (29) 

Notice that, since the filter strength depends only on this does not destruct the rotational invariance 
of the scheme Furthermore, as the order of the spherical harmonics, increases, the effect of filtering 
automatically decreases, so that the convergence of the scheme for A/" ^ 00 is retained. More specifically, 
for a filter of order p, we expect a convergence order of ~ p — 1, as suggested by Vandeven's theorem for 
Fourier expansion |46j, see also [23 . 

In our analysis we have considered two second-order and two fourth-order filters. The first one is the 
classical second-order Lanczos filteiEl 

sin 7^ 

CrLanczos(^) = , (30) 

V 

while the second and third choices are given by the ErfcLog filter [6^, 



of order p = 2,4. Finally, the fourth filter is the fourth-order spherical-spline filter 

CrsSpline(^) = ' ^^^^ 



We point out that, with our definition, the spherical-spline filter (32) is very similar to the one used in [29 



but is not exactly equivalent. The reason for using a slightly different filter form is that the filter of [29] is 



not compatible with the form of Eq. (28) as it cannot be written in terms of a function <j(-) of i/{N + 1). 

In addition, we have considered only even-order filters since the truncation error associated with these 
filters can be interpreted as a numerical viscosity of order higher than two [32], while for odd-order filters 
the leading truncation error is of the dispersion type [32 . Moreover, we also do not consider filters of orders 
higher than 4. This is because, as we will see later, the fourth-order filters are already too weak to completely 
remove oscillations, suggesting that even higher order will be even less efficient. 



^ Here we ignore the technical requirement for Vandeven's theorem that cr^^^(l) = for = 0, 1, 2, . . . p — 1, which is not 
satisfied by our filters (nor by the one proposed by 29 ). This is a condition that does not infiuence the formal accuracy of the 
filtered expansion with respect to the unfiltered truncated expansion, but it is mainly needed to prove the convergence of the 
filtered expansion in the case in which the unfiltered expansion is not converging point-wise (for instance due to the presence 
of discontinuities) [22] [45] , 

^This is a consequence of the classical addition theorem for spherical harmonics (e.g., lOl). 

^Note that the Lanczos filter is usually defined as a{r]) = sinTTTj/TTT] to have a first-order zero at 77 = 1, as discussed 
in footnote [T] Our modified Lanczos filter yields a more uniform damping of high-order modes and works very well in our 
experiments. 
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In our scheme, we filter the solution after each sub-step of the time integrator. This yields the following 
scheme 

= + e\ + S^bF\+„2. (33) 

F\+y2 = -^^bF^., (34) 
— F^f^ 

= A [i^/c+l/2]+e /e+i/2+5' bF (35) 



At 



?A_ _ TpB 



k+1 



^^bF^**, (36) 



where is a diagonal matrix representing the filtering operation. 

We should remark that both our scheme and the one by cannot be interpreted as a continuum 



problem, in the sense that the equations (33) — (36) do not, in general, represent a discretization of any system 
of partial differential equations. The main reason is that is not idempotent, i.e., ^^c^^b ^^b^ so 



that the filtering operations in Eqs. (34) and (36) do not have a well-defined behavior in the limit At 0. 
In the case in which ^^b is idempotent, the scheme has indeed a continuum limit, but it can be easily 
demonstrated that the FP^ method is just the Pm method for some M < 

This problem can be solved by making the strength parameter s (and ^^b with it) depend on the 
timestep. In order to see that, we consider the behavior of our scheme for a given mode u = F^^ with 
I ^ {). Let q = cr(l/{N + 1)), where a is a filter function. Then the effect of filtering on u in each of the two 



filtering steps in, e.^., Eq. (36) is simply: 
This can be rewritten as 

At/2 =Mj2^'' 
If we let s = (3 At, then, in the limit of At 0, we obtain 



^ = (3logqu. (37) 



In other words, we can interpret filtering as a first-order, operator split, discretization of the system of 
equations 

TP A 

T^'^'b ^ = + S^'bF'' + I3L^bF'' , (38) 



dF^ I r'^AjF'' 

dt dx^ 



where L^b is a diagonal matrix with coefficients log<j(£/(7V + 1)). This is the desired continuum limit. The 
physical interpretation is that filtering is equivalent to a forward-peaked scattering operator (notice that 
cr(0) = 1). Finally, we note that we can estimate the filter effective opacity by looking at the dissipation 
rate for the highest-order multiple of the expansion as 

(Jeff = -/3 log a(7V/(7V + l)). 

4. Tests 

In this Section, we present some tests of the numerical schemes described above as implemented in our 
Charon code. Charon uses 3D Cartesian coordinates in space and is currently parallelized employing hybrid 
OpenMP/MPI parallelization using the domain decomposition method. It uses the open-source Cactus 
Computational Toolkit [21^^, which provides MPI parallelization, input/output, and restart capabilities. 



-•^^The reason is that the only idempotent filter is the cut-off filter, that is, the filter that simply sets to zero all the modes 
with i > M for some M, while leaving unaffected all the modes with i < M. 
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4.1. ID diffusion of a step function 

In this first test, we primarily focus on verifying the abihty of our scheme to handle diffusion of radiation 
when the opacity is high and the me an- free- path is small compared to the grid spacing. In this limit, the 
continuous hyperbolic transport equation displays parabolic character to leading order [31]. Despite this, 
there is no guarantee that a numerical scheme for solving the hyperbolic system will be AP, that is, will 
reproduce a valid discretization of the asymptotic limit of the continuous equations (cf. the discussion in 



Section 3.4). 



Consider therefore the following initial conditions for the radiation energy density E = J Idftdu^ in a 
non-moving infinite medium with a constant (isotropic and elastic) scattering opacity: 

E{z, t) = R{z + l/2)H(l/2 - z) , (39) 

where H(-) is the Heaviside step function. If the scattering opacity is high, the solution of the transport 
equation is well approximated by the solution of the diffusion equation. The corresponding diffusion equation 
has the following analytic solution to problem (|39|) 




(40) 



where Erf(-) is the error function and r = 3/^s/c is the diffusion timescale, where hCg is the total scattering 
opacity, which we set = 10^ [e.^.,|3r. 

We employ five different schemes for this test: the step scheme (i.e., a DG scheme with step-limiter, or, 
equivalently a first-order FV scheme), two linear DG methods employing minmod and minmod2 limiters 
and two finite- volume (FV) methods also employing minmod and minmod2 limiters. In our implementation, 
the finite-volume scheme is obtained from the linear DG scheme by simply replacing the linear DG slope 
with the one obtained from the reconstruction procedure. 

In all of our runs, we use the Pi scheme because in ID there are no negative solutions and thus filtering 
is not necessary, and because the radiation is nearly isotropic in such a diffusive regime so that Pi scheme 
should be sufficiently accurate. We perform calculations using three different resolutions Az = 0.16,0.08 
and Az = 0.04, with the grid ranging from —2 to 2, and imposing periodic boundary conditions at the 
outer boundaries. We choose the CFL factor to be 0.25 and we recall that the maximum CFL factor that 
guarantees the L^-stability of our scheme is 1/3 in ID [17J. In all of our tests, the CFL factor is mainly chosen 
for convenience in order to have a sufficient number of timesteps within a given time interval. Moreover, 
in many radiation-transport calculations, in the absence of hydrodynamical equations, the truncation error 
due to the time discretization is expected to be small compared to other sources of error (e.^., the angular 
and spatial discretization). 

Figure [2] shows the radiation energy density as a function of z coordinate for the run with Az = 0.04 at 
time 100/c. The thick black line corresponds to the analytic solution, while the other lines show numerical 
results obtained with the above schemes. The line with red circles corresponds to the linear DG with 
minmod2 limiter solution. The line with orange squares represents the solution from the finite-volume 
scheme with minmod2 reconstruction. Finally, the lines with blue diamonds and green triangles show the 
results obtained with linear DG with minmod and step limiters, respectively. Note that the different symbols 
also mark the value of the numerical solution at each grid point (i.e., we show two points for each element). 

The linear DG method with minmod2 agrees well with the analytical result. This is expected since 
this scheme has the correct asymptotic limit [31 . In contrast, all other schemes overestimate the diffusion 
rate. In particular, the step scheme produces the worst results. It reaches stationarity already at time 
t/r ^ 10~^, which is much smaller than the diffusion timescale for this problem. The linear DG and finite- 
volume schemes with minmod yield identical results and for this reason we show only the results from the 
linear DG scheme. Both are only marginally better than the step algorithm. These results are in overall 
agreement with the ones reported by [ST for a very similar test. 

The FV scheme with the minmod2 reconstruction produces results that are relatively accurate compared 
to the linear DG and finite-volume schemes with the minmod limiter, even though the observed diffusion 



13 



0.4 



1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 

^— Diffusion Eq._ 
♦ Minmod I 
•— • l\/linmod2 
□ □ FV-l\/linmod2 - 
A— A Step - 

□ 

□ 


=■-□ □ □ □ □ □ rn V 

j 1 




~ 1 1 1 1 1 1 1 1 1 1 




1 1 1 1 1 1 1 1 1 ~ 



0.0 0.2 0.4 0.6 0.8 1.0 

Z 



Figure 2: Radiation energy density as a function of the z coordinate at time 100/c for the ID diffusion of the step func- 
tion ( |39| . The thick black hne represents the analytical solution of the corresponding diffusion equation, the line with red 
circles corresponds to the linear DG with minmod2 limiter solution, the line with orange squares represents the solution from 
the finite-volume scheme with minmod2 reconstruction, while the lines with blue diamonds and with green triangles show 
the results obtained with linear DG with minmod and step limiters, respectively. The symbols also mark the values of the 
numerical solution at each grid point (i.e., we show two points for each element). 



timescale (at the current resolution) is still orders of magnitude larger than the physical one. Finally, we 
point out that the results obtained from the other lower-resolution runs (not shown here) are in overall 
agreement with the ones presented here. 

4-2. ID diffusion of a sine wave 

Next, we explore convergence of the numerical solution to the asymptotic one. To this scope we consider 
the ID diffusion of a sine wave and thus adopt as initial conditions the energy density given by 

(41) 



(42) 



E{z, 0) = 3V47r sin (^y ^ + 1 , 
which has the following analytic solution in the diffusion limit 



E{z,t) = 3V47r 



l + exp(---lsm(-j 



For this test, our computational domain ranges from —3 to 3 and we use eight different resolutions ranging 
from = 0.3 to = 0.009375. The CFL factor is chosen to be 0.3 

Figure [3] shows the L^— norm of the deviation of the numerical result from the asymptotic solution as 
a function of numerical resolution at time t — 1000/c. As expected, the linear DG scheme with minmod2 
(line with red circles) exhibits approximately second-order convergence for the entire range of resolutions 
shown in the plot, while the linear DG with minmod (line with blue squares) starts converging only when 
Az ~ 10~^, afterwards it converges with order 1.28. The step scheme (line with green triangles) does not 
show any sign of convergence. These results are again consistent with what was observed in [31j. Finally, the 
finite-volume with minmod2 reconstruction (line with orange diamonds) exhibits second-order convergence 
even though this scheme is not asymptotic preserving. 

4.3. The 2D line-source problem 

As a first multidimensional problem used to benchmark different implementations of the filtered spherical 
harmonics discretization schemes we consider the so-called "line-source" problem, where we have initial 
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Figure 3: L°^— norm of the deviation of the numerical result from the asymptotic solution as a function of numerical resolution 
for the diffusion problem of a sine function. The lines with red circles and blue squares correspond to the linear DG schemes 
with minmod2 and minmod limiters, respectively, while the line with orange diamonds represents the FV scheme with minmod2 
reconstruction. Finally, the line with green triangles corresponds to the step scheme. 



conditions given b}pl 

/(x,^,z,l],t) = — (43) 

which represent an isotropic pulse of radiation with a total energy concentrated along the z-axis in 
vacuum [8 . This radiation field propagates in vacuum as it evolves in time according to the following 
analytical solution 

E[x, y) = ; ^ , (44) 

27r ct\J (?t^ — 



where r = ^J ^ is the distance from the z-axis and H(-) is the unit step function and (5(-) is the Dirac 
delta function. According to this solution, the radiation field consists of a front that forms a cylindrical 
shell that travels outwards at the speed of light, while in the interior, the radiation energy density smoothly 
decreases along the radial direction. 

We point out that, while this test is actually one-dimensional in cylindrical coordinates, it becomes 
particularly challenging for radiation-transport codes, except for Monte-Carlo codes, when solved on a two- 
dimensional Cartesian grid (as we do). In these coordinates, the radiation beam, which originates from a 
single spatial grid zone, has a very forward-peaked distribution in angle. This is a huge challenge for both 
spatial and angular discretization schemes. Moreover, such a configuration favors negative solutions in the 
P/v scheme. Indeed, the analytical P/v solution to this problem was shown to have regions with negative 
values of the energy [HI |30l |29], while the P\ solution even exhibits a negative singularity [27 . For all 
the results presented here, we use a grid with resolution Ax = Ay = 0.02 and a CFL factor of 0.0625. 
Furthermore, we choose Eq = \f^. 

Figure |4] displays the colormap of the radiation energy density in the x — y plane at i — 1/c. The 
upper left panel shows the analytic solution, while the upper right panel shows the pure Pj solution (note 
the considerable difference in scale). As expected, the Pj solution exhibits unphysical oscillations in the 



radial direction that are absent in the analytical solution to the full transport problem (44). The lower left 
panel of Fig. [4] shows instead the FPj solution with spherical- spline filter with effective opacity (Jeff = 20 (the 
dependence of the solution on the order N and on the value of the filter strength (Jeff will be discusses below) . 



-"^The initial conditions are 3D but we exploit the cylindrical symmetry to solve the problem on the {x^y) plane only. 
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Figure 4: Colormaps of the radiation energy density on the x — y plane at t = 1/c for the hne problem with different methods. 
The upper left panel shows the analytic solution, while the upper right panel shows the pure P7 solution (note the considerable 
difference in scale). The lower left panel shows the FP7 solution with spherical-spline filter with effective opacity CTeff = 20, 
while the lower right panel shows the FPj solution with the Lancsoz filter with the same effective opacity. 
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Figure 5: The radiation energy density as a function of the x coordinate dX t = 1/c for the hne problem test. The thick black 
line corresponds to the analytic solution, while the rest of the lines represent the results from P7 calculations without filter 
(blue line), with fourth-order ErfcLog filter (green line), with the second-order ErfcLog filter (red line), with the Lanczos filter 
(cyan line), and with the fourth-order spherical-spline filter (magenta line). 



In this case, the radial oscihations are significantly reduced compared to the unfiltered Pj solution, similarly 
to what was found in [29] . Finally, the lower right panel shows the FPj solution with the Lancsoz filter with 
the same effective opacity. In this case, the amount of oscillations is even smaller and we get a result that 
is closer to the analytical one. The reason seems to be that the Lanczos filter, being a second-order filter, 
is more effective in reducing the appearance of oscillations. The solution obtained with the spherical-spline 
filter is still characterized by the presence of a ring structure that resembles the more oscillatory (unfiltered) 
P/v solution. This structure does not disappear even for the values of the filter strength as high as deff = 10^, 
suggesting that this is a result of a shortcoming of this particular filter. We point out that we have repeated 
these runs with the second-order (fourth-order) ErfcLog filter and we obtain a result similar to the one with 
the second-order Lanczos (fourth-order ErfcLog) filter, suggesting that the order of the filter plays the most 
important role, at least for this test. 

A more quantitative measure of this test is shown in Fig. |5) where we plot a ID cut of the radiation 
energy density as a function of the x-coordinate at t = 1/c. The thick black line again corresponds to 
the analytic solution, while the rest of the lines represent the results from P7 calculations without filter 
(blue line), with the fourth-order ErfcLog filter (green line), with the second-order ErfcLog filter (red line), 
with the Lanczos filter (cyan line), and with the spherical-spline filter (violet line). All of these runs with 
filters are performed using a filter strength of (Jeff = 20. We can easily notice again the presence of large 
oscillations in the unfiltered P7 solution (which are larger than the scale of the plot). The fourth-order 
spherical-spline and ErfcLog filters are able to suppress most of the oscillations and remove the negative 
values. Nevertheless, both solutions are still affected by the oscillations (although to a much smaller extent 
compared to the P7 solution). The second-order Lacnzos and ErfcLog filters, on the other hand, are able to 
remove most of the oscillations and give the best numerical solutions, overall. 

The left panel of Fig. [6) which reports the radiation energy density as a function of the x coordinate 
at t = 1/c, highlights how the quality of the solution varies with the (Lanczos) filter effective opacity. The 
black line again corresponds to the analytic solution, while the rest of the lines represent the PP7 solution 
with the Lanczos filter of varying strength (Jeff. In the case of weak filters (e.^., deff = 1 or (Jeff = 5), there 
are significant oscillations, whose amplitude is significantly reduced as we increase the filter strength. For 
example, for <Jeff = 20, there are tiny oscillations, while for (Jeff = 50 there are no noticeable oscillations. 
However, the quality of the solution with deff = 50 is actually worse than the one with (Jeff = 20 (e.^., 
the radiation wavefront lags significantly behind the real solution) as a result of the large smearing of the 
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Figure 6: Left panel: radiation energy density as a function of the x coordinate at t = 1/c obtained from the FPr solution with 
the second-order Lanczos filter of varying strength (Jeff. The thick black line corresponds to the analytic solution, while the 
rest of the lines represent the numerical solutions corresponding to different (7^^. Right panel: same as in the left panel but 
for the FP7 solution with the fourth-order spherical-spline filter of varying strength (Jeff. 



radiation beam produced by the excessive filtering (this is even more evident for (Jeff = 100 or (Jeff = 1000). 
Therefore, the filter strength needs to be chosen large enough to damp oscillations and small enough to 
avoid excessive smearing of the solution. We have repeated these runs with the second-order ErfcLog filter 
and we again get similar results. 

The left panel of Fig. [6] should be contrasted with the right one, where we show a study of the effect of 
different filter opacities for the same problem, but using the spherical-spline filter. The first thing to notice 
is that the spherical-spline filter is never able to completely remove the "ring" structure in the solution, even 
when the filter strength is so strong that the result resembles the solution of the diffusion equation for this 
problem. Secondly, the dependence of the filter behavior on the filter strength does not seem to be easily 
predictable: at first, as we increase filter strength, negative solution disappear (for <Jeff ^ 100), then they 
reappear for higher values of cTeff around 1000. We have repeated these runs with the fourth-order ErfcLog 
filter and we obtained similar results. 

Finally, Fig. [7| shows the FPn solutions for = 3, 5, 7, 9 and 11 with the Lanczos filter with (Jeff = 20, 
and can be used to study the convergence of the FP^ approximation to the analytic solution. In this plot, 
we can distinguish three different types of errors: (1) An error in the position of the radiation front, which 
is particularly evident for small is mainly related to the fact that the propagation speed of radiation 
is smaller than c for low N. (2) An error in the profile of the radiation energy density behind the front, 
which is again particularly pronounced for small A^, and is due to the fact that high angular resolution is 
needed to properly describe the very forward-peaked angular distribution of radiation. (3) A relatively large 
spreading of the radiation beam in space compared to the x = ct — singularity in the analytical 

solution (44). This is an artefact of spatial discretization and mainly stems from the fact that the radiation 
beam originates from one spatial element and results in the presence of a "precursor" in the radiation front 
for high-order FPn {e.g., FPn) solutions, where the spatial discretization scheme propagates the radiation 
front superluminally, despite the fact that the characteristic speeds of the FPn system are always smaller 
than c. As we can see in the figure, the FPn solution nevertheless approaches the analytical one as we 
increase the order N and, in particular, the errors associated with the angular discretization decrease to the 
point that the FPn solution yields only a small improvement with respect to the FPq one. In particular, a 
large contribution to the error in the FPn solution comes from the presence of the superluminal precursor 
discussed above. Since this precursor can only be attributed to the spatial discretization error, we can 
conclude that, at this particular resolution and order, the spatial discretization error is already comparable 
with the angular discretization ones. We have repeated these runs with the second-order ErfcLog filter and 
at half of the resolution and we again get similar results. These results show that our filtering strategy is 
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Figure 7: Radiation energy density as a function of the x coordinate at t = 1/c obtained from the FP^ solution with the 
second-order Lanczos filter with (Jeff = 20 for different values of order A^. The thick black line represents that analytical 
solution, while the rest of the lines represent numerical solutions for different A^. Clearly, the FPjsf solution approaches the 
analytical one as we increase the order A^. 



able to recover the convergence of the Pn scheme for this particular case, and that the second-order filters, 
unlike the fourth-order ones, do not require a delicate fine-tuning of the effective scattering opacity. In 
particular, (Jeff can be chosen on the basis of the physics and geometry of the problem, in a way that is 
independent of the order of the employed P/v scheme. 

Overall, the results of this test confirm that the filtering approach is effective and robust in suppressing 
unphysical oscillations in the P/v solution, even with moderately low-order N . Moreover, we find that the 
second-order Lancsoz and ErfcLog filters produce significantly better numerical results compared to the 
fourth-order spherical- spline and ErfcLog filters. 



4.4. A lattice problem 

Next, we consider another 2D problem consisting of a chessboard of highly scattering and highly absorbing 
square regions located around a central emitting square region. Although this geometry is not expected to be 
present in the astrophysical scenarios we are most interested in, it nevertheless represents a very demanding 
test of the capabilities of the different numerical schemes in complicated geometries. 

In our calculation we use a setup illustrated in the upper left panel of Fig. [8) which consists of a central 
emitting square (shown in white) and 11 absorbing squares (shown in black) with a constant absorption 
opacity Ka = 10 surrounding the central emitting square. The space between the squares (shown in gray) 
and the central emitting region have a small uniform scattering opacity of = 1. Each square has a linear 
size of 1. The size of the computational domain is 7 along both axes. We choose a spatial resolution of 
Ax = = 0.035 and the CFL factor was set to ^ 0.14, with outgoing boundary conditions imposed at the 
outer boundary. 

The remaining panels in Fig. [S] show the log^g radiation energy density produced by different 

schemes at a time t = 3.2/c, which roughly corresponds to the moment when the radiation reaches the 
outer boundary of the computational domain. The upper right panel corresponds to the P7 solution. Not 
surprisingly, this solution has regions of negative energy density (shown in white), although the negative 
values reach at most a relatively small magnitude of ^ 10~^. All other solutions are computed with the 
Lanczos filter with effective opacity (Jeff = 5. We find that deff ^ 5 is necessary to avoid the appearance of 
negative solutions. The middle-left panel represents the PPi solution, which does not have negative regions, 
but where the radiation wavefront has reached only half of the computational domain. This is again d ue t o 
the fact that the N = 1 wave travels at a slower velocity of c/a/3 [37 (see the discussion in Section 3.2). 
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Figure 8: Upper-left panel: illustration of the setup of the 2D lattice problem. The rest of the panels: colormap of the log^o of 
the radiation energy density at time t = 3.2/ c as obtained with the P7 scheme (upper-right panel) and FPjy schemes with the 
Lanczos filter with opacity 5 for different order N (middle and bottom panels). The time t = 3.2 /c corresponds to the moment 
when the radiation front first reaches the outer boundary of the computational domain. 
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Figure 9: Colormaps of the log^^Q of the radiation energy density at time t = 16/c obtained with the P7 scheme (upper- left 
panel) and FPi (lower-left panel), FP7 (upper-right panel), and FP3 (lower-right panel) schemes with the Lanczos filter with 
opacity 5. The time t = 16/c corresponds to the time by which the radiation field reaches the stationary state. Since the Pn 
scheme is less likely to exhibit negative solutions in the stationary state, we do not observe such solutions in our numerical 
result (upper-left panel), similarly to the filtered Pn solutions. 



The middle- right panel shows the FP3 solution, which also does not have regions with negative energy 
density and where the front has travelled enough to cover ~ 95% of the computational domain, as a result 
of larger propagation speeds with higher TV (of course, the velocity is always bounded by the speed of light). 
The FP5 solution shown in the lower-left panel is very similar to the FP3 case, with the only noticeable 
difference being the slightly faster propagation velocity in the PP5 case. Finally, the lower-right panel shows 
the PP7 solution, which looks almost indistinguishable from the PP5 one, suggesting that the PP5 scheme 
yields sufficiently accurate results for this problem. 

We complete the analysis of this test by showing in Fig. |9] equivalent snapshots of the energy density at 
a later time of t = 16/c, when the radiation field has reached a stationary state. In these conditions, the Pat 
solution is not expected to have any negative values [8 . This is indeed confirmed by the upper left panel 
of Fig. |9j which shows the P7 solution without negative regions. The remaining panels in Fig. [9] report the 
PPi, PP3, and PP7 solutions, respectively. Note that the PP3, and PP7 solutions are very similar to the P7 
one, underlining that the filter we use does not compromise the accuracy of the solution. However, the PPi 
solution appears to be significantly different from the PP3 and PP7 solutions, implying that = 1 is not a 
sufficiently accurate approximation for this problem. 
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Figure 10: Radiation energy density as a function of the radial coordinate for the homogeneous sphere problem. The thick 
black line shows the analytical solution, while the rest of the lines show the FP^ solution with the Lanczos filter with (Jeff = 1 
but with different values of order N. Clearly, the numerical result approaches the analytical solution as we increase N. 



4-5. 3D Homogeneous sphere 

Finally, we consider the 3D homogeneous sphere problem, which is frequently employed to test radiation 
transport codes [43l |4ll [2] . This problem consists of a static homogeneous and isothermal sphere of radius 
R that radiates in vacuum. Inside the sphere, the radiation interacts with the background matter only via 
isotropic absorption and thermal emission. Despite the rather simple setup, the sharp discontinuity at the 
surface of the sphere is a model for astrophysical phenomena with rapidly varying opacity. This represent a 
major challenge for finite-difference methods (although, it is less challenging for Monte Carlo methods; see, 
e.g., 12]). 

We assume that the sphere of radius R has a constant absorption opacity hia and emissivity B in the 
interior, while in the ambient vacuum at r > i?, we have Hia = B = 0. For this problem, the transport 
equation can be solved analytically and has solution \43\ 



B 



[1 - exp(-Kas(r, ii))] 



(45) 



r=0 



where r = y^x^ + + z^, fi = cos and 



s(r, /i) 



rji^Rgir^ji) if r < R^ — 1 < < 1 ; 



2Rg{r,fi) 




if r>R, 
otherwise. 



1 - {R/rf < /i < 1 , 



(46) 



and 



(47) 



Note that this solution depends only on three parameters: R^ and 5, where the latter acts as a scale 
factor for the solution. 

We perform simulations in full 3D with Cartesian coordinates and use the following computational 
setup. We set = 1 and cover the interior of the sphere with 40 elements in diameter along each coordinate 
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Figure 11: L^— norm of the deviation of the FP/v solution from the analytic result as a function of the order for the 
homogeneous sphere problem. The black dots show the error as computed in a sphere of radius R = 4.5, while the black line 
shows the fitted convergence rate. 



direction, with the outer boundary being located at 5 R. The absorption opacity is chosen to be tZa = 10 
and the CFL factor is set to be ^ 0.12. 

It is useful to remark that, although the matter distribution is spherically symmetric, this is a genuinely 
3D test due to the Cartesian geometry of our spatial grid. Indeed, it leads to the propagation of radiation 
from one spatial zone to another not only in the radial direction, but also in the angular directions, and the 
degree of sphericity of the numerical solution can be taken as a measure of the accuracy. 

Figure 10 shows the radiation energy density along the diagonal direction for the analytical solution 
and the FPn solutions of different orders ranging from 1 to 11. The results are shown for the time when 
the radiation field has reached a stationary state{^ These runs are performed with the Lanczos filter of 
(Jeff = 1. Interestingly, all of the FP^ solutions produce the correct result in the interior of the sphere. This 
is not surprising since the radiation inside the sphere is nearly isotropic and the low-order FPn solutions 
are already accurate enough. Outside the sphere, radiation streams freely outwards with a highly forward- 
peaked distribution in angle, which is a challenge for low-order FPn schemes. Indeed, the FPi result 
deviates significantly from the analytical solution in that region. However, the solution clearly becomes 
more accurate everywhere in the computational domain as we increase the order of the scheme. Figure [TT] 
shows the L^— norm of the deviation of the FPn solution from the analytic result in a sphere of radius 
R = 4.5j^ As we can see from the plot, the FPn scheme starts to converge already for ^ 3, an order 
with only 4^ = 16 angular degrees of freedom. The convergence order is o:^ 1.16, which is consistent with 
what expected from the theory of spectral filtering. 

Figure 12 shows the radiation energy density along the diagonal direction for the analytical solution 
and different numerical solutions at stationarity state. These are the unfiltered P7 solution, the three 
FPy solutions computed using respectively the second-order Lanczos, the fourth-order spherical-spline and 
ErfcLog filters with (Jeff = 1. Also shown are the two FP7 solutions computed using the Lanczos filter with 
(Jeff = 0.1 and (Jeff = 10. Although the P7 solution does not exhibit any negative solutions, it shows large 
oscillations in the free streaming region. The FP7 solution with the Lanczos filter with (Jeff = 0.1 also yields 
a somewhat oscillatory solution, suggesting that the filter effective opacity is too low for this problem. As 



■"^^Note that stationarity is reached at different times depending on the different schemes used. For this reason and given the 
high computational costs, we did not evolve all the models up to the same time. Instead we report the solution as obtained as 
soon as stationarity is reached. In all cases the computations are performed up to at least t = 20/c. 

-"^"^We compute the error inside R = 4.5 instead of = 5 in order to exclude effects due to boundary conditions. We also 
normalize the L^— norm by dividing it by |7rfl3. 



23 




Figure 12: Radiation energy density as a function of the r coordinate for the homogeneous sphere problem. The thick black 
line shows the analytical solution, the blue line corresponds to the unfiltered P7 solution, while the rest of the lines represent 
the FP7 solutions obtained with different filters and different values of (Jeff • 



(Jeff is increased, the spurious oscillations disappear and all the filters that we have tried yield solutions of 
very similar quality for deff = 1. Finally, the FP7 solution with the Lanczos filter with (Jeff = 10 is similar 
in quality to the FPs solution with the same filter but with (Jeff = 1. This is due to the excessive damping 
of the high-order multiples of the solution by the filtering procedure. 

Finally, Fig. [Ts] shows the colormaps of the log^g radiation energy density from the FPi (left panel) 

and FPii (right panel) solutions with the Lanczos filter with (Jeff = 1 at t = 3.75/c, which corresponds to 
the time when the radiation front almost reaches the outer boundary. Clearly, and as observed also in 
the previous tests, the the FPi radiation front lags significantly behind the FPu one because of its slower 
propagation speed (cf. the discussion in Section 3.6). Moreover, we can also see that both solutions maintain 
a high level of spherical symmetry despite the Cartesian geometry of the spatial grid. We obtain similar 
results in tests where the sphere was covered by 20, 10 and 5 elements. 



5. Conclusion 

We have presented an extension of the filtered spherical harmonics method by McClarren and Hauck 
[29] . the FP/v scheme, to three dimensions. We have developed the new 3D/multigroup radiation trans- 
port code Charon, built within the Cactus Computational Toolkit [2TJ[T]. Charon uses an asymptotic- 
preserving linear discontinuous Galerkin discretization scheme in space [31] and a semi-implicit time inte- 
gration scheme [28] (cf. Section |3|. 

Our filtering scheme differs from the one presented by [29] in one important aspect: we reformulate the 
filtering procedure so that it acquires a well-defined continuum limit. In particular, we have shown that in 
the limit where the spatial and time steps are reduced to zero, our filtering scheme can be interpreted as the 
addition of a forward-peaked artificial scattering term to the P/v equations. The filtering procedure is also 
constructed in such a way as to retain the convergence of the FPn solution to the solution of the transport 
equation as TV ^ 00. 

We have tested our scheme against a few challenging benchmark problems for radiation transport using 
four different filtering kernels: the fourth order spherical- spline filter, which is similar in spirit to the filter 
used by [29], the fourth-order and second-order ErfcLog filters [6], and the classical second-order Lanczos 
filter [6]. Our findings indicate that the FP/v scheme behaves well also in the three-dimensional case. In 
addition, we have shown that the second-order filters are more robust and accurate and require somewhat 
less tuning of the filter strength when compared to the fourth-order spherical-spline and the ErfcLog filters. 
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(a) FPi solution (b) FPn solution 

Figure 13: Colormaps of the radiation energy density at time t = 3.75/ c obtained with the FPi and FP7 schemes using the 
Lanczos filter with opacity (Jeff = 1- The time t = 3.75/c corresponds to the moment when the radiation front almost reaches 
the outer boundary of the computational domain. 

Since the order of a filter is one of its most important properties, this result is likely to apply also to several 
other second- and fourth-order filters. 

In future work, we plan to extend our numerical algorithms to include velocity dependence, the coupling 
to hydrodynamics, and, eventually, general relativity. 
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Appendix A. Real spherical harmonics 

This Appendix is dedicated to the derivation of the real spherical harmonics, whose implementation in 
Charon has been particularly advantageous. We start by recalling that the spherical harmonics, are the 
eigenfunctions of the Laplace-Beltrami operator. A, on the unit 2-sphere: 

Spherical harmonics are usually written, in complex form, as 

YP {if, 0)=e''^'^PP {cos 0), 
where —i<m<i^ Pp are the associated Legendre function, see e.g. [7], 

Pr{x) = {l-x'r/'CZVi^'{x), m>0, 
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and are the Gegenbauer polynomials of index a and degree n. We, also, use the standard convention 
that 

p,-(x) = (-ir|^prw- (A.i) 

The associated Legendre functions of index m > are orthogonal in [—1, 1] with unit weight 

•'frwpywci:.- p,y,;r_",„)/ «.. 

while the associated Legendre functions of degree i and index m^m' > are orthogonal in [—1,1] with 
weight (1 -x^)-^ 



x: 



Prix) Prix) 



dx 



(£ + m)! 



1 — x'^ m{i — m)\ 

The spherical harmonics with index m, > are then orthogonal with unit weight on the sphere: 

/ Yp{ip,0)Y;r'\ip,0)dn= e'^'^-'^'^^dip PP {cos 0)Pf {cos 0)smOdO 
Jsi Jo Jo 
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i{m — m')Lp 



dip / Prix)Pr ix)dx 



_ 2ie + m)\ ^ 

~ (2^+l)(£-m)! " 

We can then redefine the spherical harmonics as 



(A.2) 
(A.3) 
(A.4) 



Yri^,0) 



' (2£ + l) ie-m)\ 
Att (£ + m)! 



e™*^ P7*(cos 9) = NP e™^ P7'(cos 9) . 



It is easy to see, that for all m,m', thanks to the normalization and the convention (A.I), we have 

0) ((^, 0) dQ = 5mm' hi' . 



Si 

We can construct a real basis from the spherical harmonics by defining 

^^(r^- + (-i)-r^-"^), ifm>o, 

Yim = \ y^, ifm = 0, 

^^(r,— -(-irr-) ifm<o. 

Using again the fact that, N'^"^^ P'^"^^ = (-1)1^1 Arl^'pj"^' , we obtain 

Yim{^,0) = V2NP cos(m(/:))Pf (cos6>), 
Yem{^,0) = V^TV]"^! sin(|m|(^) P'"^' (cos^). 



m > 0, 
m < 0, 



which is the wanted expression for the real spherical harmonics. 

We can construct a Gaussian quadrature associated with the spherical harmonics basis as a direct product 
of a uniform quadrature in the (f direction: 



M^ + 1/2' 



M^ + 1 



m. 
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which is accurate to order 4:M^ + 2 [7 and a Gauss-Legendre grid in the 6 direction [T2] : 

2 

/i = cos6>, we = — lie = {ZIPmo+iM]} p < £ < Me , 

where Pm is the Legendre polynomial of degree M, i.c.Pm = Pm ^[p] denotes the set of the roots of p. 
This quadrature formula is then accurate up to order 2Mq-\-1. This means that if we want to obtain an exact 
representation of the scalar product of spherical harmonics up to order Nq we have to choose Mq = Nq. 
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